home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_4.4 / XPM / PROC_POL.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  8KB  |  280 lines
  1. /* 
  2.  * proc_pol.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * PROC(ess)_POLY(gon)
  11.  *
  12.  * polygon analysis:
  13.  * -----------------
  14.  *   reconstruct polygon;
  15.  *   determine area by pixel counting
  16.  *            (and, in poly_moments(), by evaluating zero order moment);
  17.  *
  18.  *   evaluate curvature energy based on Zahn chain code
  19.  *
  20.  *   compute centroid, central moments and principal axes;
  21.  *
  22.  *   find 2**n new vertices equi-spaced along contour
  23.  *   determine radius vector relative to centroid; shift to zero mean;
  24.  *   reconstruct angular bend function with respect to new set of vertices;
  25.  *
  26.  *   compute power_spectrum and autocorrelation by AP FFT,
  27.  *           for angular bend function (-->tan mode) 
  28.  *           for radial function (-->rad mode)
  29.  *           for both (-->all)
  30.  *   save to file (writes file for plotting by routine acmp);
  31.  *
  32.  *   evaluate Zahn-Roskies Fourier descriptors
  33.  *   write them to file (writes file for plotting by modeplt, multiplt);
  34.  *
  35.  */
  36.  
  37. #include <stdio.h>
  38. #include <stdlib.h>
  39. #include <math.h>
  40. #include "ip.h"
  41.  
  42. #define    ON_MODE        1
  43. #define    OFF_MODE    0
  44. #define    C_HEIGHT    10
  45. #define    C_WIDTH        10
  46. #define    C_COLOR        255
  47.  
  48.  
  49. #define N2_LIMIT    64L
  50. #define    MAX_ZRFD    15L
  51. #define N_MOMENTS    15
  52. #define N_MODE_PARMS    2
  53.  
  54. #define ENLARGE        1.0
  55. #define NOLOOP         0
  56. #define MAX_ROW     2048
  57. #define LOOP_COMPLETE    8.0
  58.  
  59. #define SQ2        1.414213562
  60. #define ZERO        0.0
  61. #define FABS(a)        ((a) > ZERO ? (a) : - (a))
  62.  
  63. #define ON        1
  64. #define OFF        0
  65. #define SHIFT_CENTROID    ON
  66. #define    SAVE_SHAPE    OFF
  67.  
  68. #define    EVAL_CURV_EN    OFF        /* evaluate curv. energy */
  69. #define    POLY_FIT    ON             /* invoke Wall-Danielsson fit */
  70.  
  71. #define DISPLAY_DATA
  72.  
  73.  
  74. /* globals */
  75. extern int page;
  76.  
  77.  
  78. struct polygon *
  79. select_poly (struct polygon *poly_head, Image * imgIO, int value)
  80. {
  81.   struct polygon *current_poly;
  82.  
  83.   int c;
  84.  
  85.   int rad_mode = OFF, tan_mode = OFF;
  86.  
  87.   float *s, *s_ap, *s_n2;
  88.   float *moments;
  89.   float *phi_k;
  90.   float *r_n2;
  91.  
  92.   float *n2_phi_star;
  93.   float *acf, *p_spec;
  94.   int *mode_parms;
  95.   float *a_n = NULL, *b_n = NULL;
  96.  
  97.   long nv, nv2;
  98.   int n_mom = N_MOMENTS;
  99.   int n_mp = N_MODE_PARMS;
  100.  
  101.   unsigned int pixel_count;
  102.  
  103.   double av_dirn = 0.0;
  104.   double c_len, ap_c_len, res_c_len;
  105.   double tot_phi;
  106.   double curv_en, e_bend;
  107.   double p2a_ratio;
  108.  
  109.   struct Bdy bd, *bdp = &bd;
  110.   struct spoint *hullctr;
  111.   struct spoint cursor, *pc = &cursor;
  112.   struct spoint *v, *v_ap, *v_n2, *v_n2h;
  113.   struct spoint *pvc;
  114.   int hull_area;
  115.   char in_buf[IN_BUF_LEN];
  116.  
  117. /*
  118.  * cycle through list of boundary polygons
  119.  */
  120.   printf ("...enter n to cycle, y to select poly...\n");
  121.   current_poly = poly_head;
  122.   for (;;) {
  123.     pc->x = current_poly->first_x;
  124.     pc->y = current_poly->first_y;
  125.     printf ("poly starting at (%d, %d) select?(y/n) ", pc->x, pc->y);
  126.     if ((c = readlin (in_buf)) == 'y')
  127.       break;
  128.     if ((current_poly = current_poly->next_poly) == NULL) {
  129.       printf ("At end of polygon list!!  Return to first polygon? (y/n)");
  130.       if ((c = readlin (in_buf)) == 'n')
  131.         return (current_poly);
  132.       else
  133.         current_poly = poly_head;
  134.     }
  135.   }
  136.  
  137.   nv = current_poly->poly_points;
  138.   nv2 = 2L;
  139.   while (nv2 < nv)
  140.     nv2 *= 2L;
  141.   if (nv2 >= N2_LIMIT)
  142.     nv2 = N2_LIMIT;
  143.  
  144.  
  145. /*
  146.  * allocate memory 
  147.  */
  148.   if ((v = (struct spoint *) malloc ((nv + 1) * sizeof (struct spoint))) == NULL)
  149.       exitmess ("\n...mem allocation for v failed\n", 1);
  150.   if ((v_ap = (struct spoint *) malloc ((nv + 1) * sizeof (struct spoint))) == NULL)
  151.       exitmess ("\n...mem allocation for v_ap failed\n", 1);
  152.   if ((phi_k = (float *) malloc ((nv + 1) * sizeof (float))) == NULL)
  153.       exitmess ("\n...mem allocation for phi_k failed\n", 1);
  154.   if ((s = (float *) malloc ((nv + 1) * sizeof (float))) == NULL)
  155.       exitmess ("\n...mem allocation for s failed\n", 1);
  156.   if ((s_ap = (float *) malloc ((nv + 1) * sizeof (float))) == NULL)
  157.       exitmess ("\n...mem allocation for s_ap failed\n", 1);
  158.  
  159.  
  160.   if ((moments = (float *) malloc (n_mom * sizeof (float))) == NULL)
  161.       exitmess ("\n...mem allocation for moments failed\n", 1);
  162.  
  163.  
  164.   if ((v_n2 = (struct spoint *) malloc ((nv2 + 1) * sizeof (struct spoint))) == NULL)
  165.       exitmess ("\n...mem allocation for v_n2 failed\n", 1);
  166.   if ((v_n2h = (struct spoint *) malloc ((nv2 / 2) * sizeof (struct spoint))) == NULL)
  167.       exitmess ("\n...mem allocation for v_n2h failed\n", 1);
  168.   if ((r_n2 = (float *) malloc (nv2 * sizeof (float))) == NULL)
  169.       exitmess ("\n...mem allocation for r_n2 failed\n", 1);
  170.   if ((s_n2 = (float *) malloc ((nv2 + 1) * sizeof (float))) == NULL)
  171.       exitmess ("\n...mem allocation for s_n2 failed\n", 1);
  172.  
  173.   if ((n2_phi_star = (float *) malloc (nv2 * sizeof (float))) == NULL)
  174.       exitmess ("\n...mem allocation for n2_phi_star failed\n", 1);
  175.  
  176.   if ((acf = (float *) malloc (nv2 * sizeof (float))) == NULL)
  177.       exitmess ("\n...mem allocation for acf failed\n", 1);
  178.   if ((p_spec = (float *) malloc (nv2 * sizeof (float))) == NULL)
  179.       exitmess ("\n...mem allocation for p_spec failed\n", 1);
  180.   if ((mode_parms = (int *) malloc (n_mp * sizeof (int))) == NULL)
  181.       exitmess ("\n...mem allocation for mode_parms failed\n", 1);
  182.  
  183. /*
  184.  * evaluate arc_length and Zahn angular bend function
  185.  */
  186.   c_len = arc_length (current_poly->d_l_ptr, s, nv);
  187.  
  188.   tot_phi = angular_bend (current_poly->d_phi_ptr, phi_k, nv);
  189.  
  190. /*
  191.  * reconstruct polygon vertices
  192.  */
  193.   pixel_count = reconstruct_poly (current_poly->d_phi_ptr,
  194.                                   current_poly->d_l_ptr, nv, v, pc,
  195.                                   current_poly->first_edge_out,
  196.                                   current_poly->total_delta_phi,
  197.                                   imgIO, value);
  198.  
  199.   printf ("\n...object area (pixel count): %u", pixel_count);
  200.  
  201.  
  202. /*
  203.  * evalate curvature energy from chain code representation
  204.  * note:
  205.  *   to obtain robust estimates of curv_en, evaluation should really
  206.  *   not be attempted before smoothing, e. g. via polygonal approx.,
  207.  *   see below;
  208.  */
  209.   if (EVAL_CURV_EN == ON) {
  210.     curv_en = curvature_energy (current_poly->d_phi_ptr,
  211.                                 current_poly->d_l_ptr, nv);
  212.     printf ("\n...chain-code derived curv. energy: %f\n", curv_en);
  213.   }
  214.  
  215. /*
  216.  * Wall-Danielsson polygonal approximation
  217.  */
  218.   res_c_len = p2a_ratio = e_bend = -1.0;
  219.   if (POLY_FIT == ON) {
  220.     fill_bdy_structure (bdp, v, nv, tot_phi);
  221.     hullctr = poly_hull (bdp, v_ap, av_dirn, &hull_area, imgIO, GRAY);
  222.     printf ("\tlocation of convex hull center:");
  223.     printf (" (%3d,%3d)\n", hullctr->x, hullctr->y);
  224.  
  225.     if ((nv - bdp->an) > 0) {
  226.       printf ("\n...approx poly size: %ld\n", nv = bdp->an);
  227.  
  228.       v_ap = (struct spoint *) realloc (v_ap, (nv + 1) * sizeof (struct spoint));
  229.  
  230.       nv2 = 2L;
  231.       while (nv2 < nv)
  232.         nv2 *= 2L;
  233.     }
  234.  
  235. /*
  236.  * evaluate polygon moments
  237.  */
  238.     pvc = poly_moments (nv, v_ap, moments, imgIO, GRAY);
  239.     printf ("\n...centroid pos: (%3d, %3d)\n", pvc->x, pvc->y);
  240.  
  241. /*
  242.  * evaluate arc_length (partial sums of edge lengths) from set of
  243.  * vertices for approximated polygon;
  244.  */
  245.     ap_c_len = vert_to_clen (v_ap, s_ap, nv);
  246.  
  247. /*
  248.  * resample new polygon at 2**n sites
  249.  */
  250.     n2_vert (v_ap, s_ap, nv, v_n2, nv2, imgIO, 192);
  251.   }
  252.   else {
  253.     pvc = poly_moments (nv, v, moments, imgIO, GRAY);
  254.     n2_vert (v, s, nv, v_n2, nv2, imgIO, 192);
  255.   }
  256.  
  257.   printf ("\n-->polygon analysis completed<--\n");
  258.  
  259. /*
  260.  * deallocate memory 
  261.  */
  262.   free (v);
  263.   free (v_ap);
  264.   free (v_n2);
  265.   free (v_n2h);
  266.   free (s);
  267.   free (s_ap);
  268.   free (s_n2);
  269.   free (phi_k);
  270.   free (moments);
  271.   free (r_n2);
  272.  
  273.   free (n2_phi_star);
  274.  
  275.   free (acf);
  276.   free (p_spec);
  277.  
  278.   return (current_poly);
  279. }
  280.